** The following do file is used to estimate the bootstrapped standard errors
** for the point estimates of the discontinuous change in cooling-driven
** electricity consumption caused by Title-24. The following steps are performed
**
** 1. From the two-year panel of daily, premise-level consumption, we generate
** a bootstrap sample by drawing 104 weeks (with replacement). This ensures that
** any correlation in consumption within, and across, households in a given week
** is maintained.
**
** 2. For the bootstrap sample, we estimate the premise-specific temperature
** response functions and the average annual cooling-driven electricity
** consumption from each premise (specified by Eq. 1 and 3).
**
** 3. From the set of premise-specific annual cooling estimates, we randomly
** draw 614 Census Block Groups (with replacement from the 614 Census Block
** Groups in our data) and estimate the discontinuous change in cooling driven
** consumption and the discontinuous change in the minimum consumption
** temperature using our primary model specified by Eq. 5.
**
** 4. We repeat steps 1-3 above 10,000 times to generate a distribution of point
** estimates for each of the coefficients in Eq. 5. The standard errors reported
** in the main RD results (Tables 2 and A9) now represent the sample standard
** deviation from the 10,000 point estimates.


set more off
set matsize 11000
clear all

cd \Data

set obs 2003
gen temp_avg = (_n + 5198)/100
mkspline temp = temp_avg, cubic knots(52 62 72)
keep temp1 temp2
mkmat temp1 temp2, matrix(temp)


mata:
void get_minindex(string scalar m, string scalar idx)
{
    real matrix A
    real scalar i
    real rowvector res
    real colvector v
    
    A = st_matrix(m)
    
    /*the min index of column i is stored in res[1, i]*/
    res = J(1, cols(A), .)
    
    /*get the min index*/
    for(i=1; i<=cols(A); i++){
        minindex(A[.,i], 1, v=., .)
        res[1, i] = v[1, 1]
    }
    
    /*save result into a Stata matrix idx*/
    st_matrix(idx, res)
}
end


** This section of the do file creates a new consumption dataset which only
** includes premises built between 1968 and 1989, have no missing values for
** the Assessor control variables that will be used in the RD models, and are
** zoned for single family residences. For this subset of the premises, a new,
** sequential identifier (id) will be created -- which ulitmately will be looped
** over when estimating the annual cooling at each premise.

use cons_data.dta, clear
sort premise_id
merge m:1 premise_id using premise_characteristics.dta

use premise_characteristics.dta, clear
keep if _merge == 3
keep if year_built >= 1968 & year_built <= 1989
keep if electric_heat == 0
keep if (zoning == "R-1" | zoning == "RD 1" | zoning == "RD 2" | zoning == "RD 3" | zoning == "RD 4" | zoning == "RD 5" | zoning == "RD-1" | zoning == "RD-2" | zoning == "RD-3" | zoning == "RD-4" | zoning == "RD-5" | zoning == "RD1" | zoning == "RD2" | zoning == "RD3" | zoning == "RD4" | zoning == "RD5" | zoning == "RD5PD")
drop if bedrooms_num == 0 | bedrooms_num == .
drop if Resdnc_sqft == 0 | Resdnc_sqft == .
sort premise_id
save Temp/independent_variables_short.dta, replace

use cons_data.dta, clear
sort premise_id
merge m:1 premise_id using independent_variables_short.dta
collapse year_built, by(premise_id)
drop year_built
gen id = _n
sort premise_id
save Temp/short_premise_ids.dta, replace

use cons_data.dta, clear
sort premise_id
merge m:1 premise_id using Temp/short_premise_ids.dta
keep if _merge == 3
save Temp/bootstrap_cons_data_short.dta, replace


** The following subsection of the do file predicts the annual cooling at each
** premise in 10,000 different bootstrap samples which are createtd by drawing
** 104 weeks worth of data with replacement.

forvalues s = 1/10000 {

use Temp/bootstrap_cons_data_short.dta, clear
bsample 104, cluster(week_id)
save Temp/pseudo_bootstrap_cons_data.dta, replace

forvalues j = 1/96 {
	use Temp/pseudo_bootstrap_cons_data.dta, clear
	keep if id > (`j' - 1) * 500 & id <= `j' * 500
	save Temp/pseudo_bootstrap_cons_data_`j'.dta, replace
}

use Temp/pseudo_bootstrap_cons_data.dta, clear
keep if id > 48000
save Temp/pseudo_bootstrap_cons_data_97.dta, replace

forvalues j = 1/96 {
	local id_min = 500 * (`j' - 1) + 1
	local id_max = 500 * `j'
	
	use Temp/pseudo_bootstrap_cons_data_`j'.dta, clear
	
	forvalues i = `id_min'/`id_max' {
		quietly: reg kwh_daily temp1 temp2 if id == `i'
		matrix define beta_`i' = e(b)
	}

	matrix define subset_`j' = (`id_min', beta_`id_min')
	local k = 500 * (`j' - 1) + 2
	
	forvalues i = `k'/`id_max' {
		matrix define subset_`j' = (subset_`j' \ `i', beta_`i')
	}
	
	clear
	set obs 2003
	gen temp_avg = (_n + 5198)/100
	mkspline temp = temp_avg, cubic knots(52 62 72)
	keep temp1 temp2
	mkmat temp1 temp2, matrix(temp)

	matrix define A = temp * subset_`j'[1..., 2..3]'
	mata:get_minindex("A", "minA")
	matrix define min_temp = (minA' + vecdiag(I(500))'*5198)/100

	clear
	svmat subset_`j'
	svmat min_temp
	rename subset_`j'1 id
	rename subset_`j'2 temp1
	rename subset_`j'3 temp2
	rename subset_`j'4 cons
	rename min_temp1 min_temp
	save Temp/subset_`j'.dta, replace
	clear
	clear matrix
}

use Temp/pseudo_bootstrap_cons_data_97.dta, clear

forvalues i = 48001/48411 {
	quietly: reg kwh_daily temp1 temp2 if id == `i'
	matrix define beta_`i' = e(b)
}

matrix define subset_97 = (48001, beta_48001)
forvalues i = 48002/48411 {
	matrix define subset_97 = (subset_97 \ `i', beta_`i')
}

clear
set obs 2003
gen temp_avg = (_n + 5198)/100
mkspline temp = temp_avg, cubic knots(52 62 72)
keep temp1 temp2
mkmat temp1 temp2, matrix(temp)

matrix define A = temp * subset_97[1..., 2..3]'
mata:get_minindex("A", "minA")
matrix define min_temp = (minA' + vecdiag(I(411))'*5198)/100

clear
svmat subset_97
svmat min_temp
rename subset_971 id
rename subset_972 temp1
rename subset_973 temp2
rename subset_974 cons
rename min_temp1 min_temp
save Temp/subset_97.dta, replace
clear
clear matrix

use Temp/subset_1.dta
forvalues j = 2/97 {
	append using Temp/subset_`j'.dta
}

save Temp/Pseudo_Sample_Estimates_`s'.dta, replace

}

********************************************************************************
** The following subsection of the code predicts the annual cooling consumption
** using the temperature response functions estimated above.

use cons_data.dta, clear
collapse (mean) temp_avg, by(year month day)
mkspline temp = temp_avg, cubic knots(52 62 72)
mkmat temp1-temp2, matrix(spline)
mkmat temp_avg, matrix(temp_avg)

matrix define ones = vecdiag(I(731))'

forvalues s = 1/10000 {

forvalues j = 1/96 {

	local id_min = 500 * (`j' - 1) + 1
	local id_max = 500 * `j'
	
	use Temp/Pseudo_Sample_Estimates_`s'.dta, clear
	keep if id >= `id_min' & id <= `id_max'
	mkspline min = min_temp, cubic knots(52 62 72)
	mkmat min1-min2, matrix(min_temp_spline)
	mkmat temp1-temp2, matrix(coef)
	
	forvalues i = `id_min'/`id_max' {
	
		local row = `i' - 500*(`j' - 1)
	
		matrix define min_temp_spline_vec = ones*min_temp_spline[`row',1..2]
		matrix define min_temp_vec = ones*min_temp_spline[`row',1]
		matrix define cooling = (spline - min_temp_spline_vec) * coef[`row',1..2]'

		svmat cooling
		svmat min_temp_vec
		svmat temp_avg
		keep if temp_avg1 > min_temp_vec1
		collapse (sum) cooling1, by(min_temp_vec1)
		mkmat min_temp_vec1-cooling1, matrix(cooling_`i')
		clear
	}

	matrix define subset_`j' = (`id_min', cooling_`id_min')
	matrix drop cooling_`id_min'
	
	local k = `id_min' + 1
	
	forvalues i = `k'/`id_max' {
		matrix define subset_`j' = (subset_`j' \ `i', cooling_`i')
		matrix drop cooling_`i'
	}
	
	clear
	svmat subset_`j'
	matrix drop subset_`j'
	rename subset_`j'1 id
	rename subset_`j'2 min_temp
	rename subset_`j'3 cooling
	save Temp/Cooling_`j'.dta, replace
	
}

use Temp/Pseudo_Sample_Estimates_`s'.dta, clear
keep if id > 48000
mkspline min = min_temp, cubic knots(52 62 72)
mkmat min1-min2, matrix(min_temp_spline)
mkmat temp1-temp2, matrix(coef)

forvalues i = 48001/48411 {

	local row = `i' - 48000

	matrix define min_temp_spline_vec = ones*min_temp_spline[`row',1..2]
	matrix define min_temp_vec = ones*min_temp_spline[`row',1]
	matrix define cooling = (spline - min_temp_spline_vec) * coef[`row',1..2]'

	clear
	svmat cooling
	svmat min_temp_vec
	svmat temp_avg
	keep if temp_avg1 > min_temp_vec1
	collapse (sum) cooling1, by(min_temp_vec1)
	mkmat min_temp_vec1-cooling1, matrix(cooling_`i')
}

matrix define subset_97 = (48001, cooling_48001)
matrix drop cooling_48001
	
forvalues i = 48002/48411 {
	matrix define subset_97 = (subset_97 \ `i', cooling_`i')
	matrix drop cooling_`i'
}
	
clear
svmat subset_97
matrix drop subset_97
rename subset_971 id
rename subset_972 min_temp
rename subset_973 cooling
save Temp/Cooling_97.dta, replace
	


use Temp/Cooling_1.dta, clear
forvalues j = 2/97 {
	append using Temp/Cooling_`j'.dta
}

sort id
save Temp/Pseudo_Sample_Cooling_`s'.dta, replace

}

********************************************************************************
** RD estimates from psuedo samples

forvalues s = 1/10000 {

use Temp/independent_variables_short.dta, clear
merge 1:1 id using Temp/Pseudo_Sample_Cooling_`s'.dta
keep if _merge == 3
replace cooling = 0 if min_temp > 72
gen cooling_kwh = cooling/2

bsample 614, cluster(block_group)

xtset comm_block

quiet: reg cooling_kwh post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=52
matrix define beta = e(b)
matrix define beta_1 = beta[1,1..3]

quiet: reg cooling_kwh post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52
matrix define beta = e(b)
matrix define beta_2 = beta[1,1..3]

quiet: reg cooling_kwh post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52
matrix define beta = e(b)
matrix define beta_3 = beta[1,1..2]

quiet: xtreg cooling_kwh post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=52, fe
matrix define beta = e(b)
matrix define beta_4 = beta[1,1..3]

quiet: xtreg cooling_kwh post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52, fe
matrix define beta = e(b)
matrix define beta_5 = beta[1,1..3]

quiet: xtreg cooling_kwh post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52, fe
matrix define beta = e(b)
matrix define beta_6 = beta[1,1..2]

quiet: reg min_temp post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72
matrix define beta = e(b)
matrix define beta_7 = beta[1,1..3]

quiet: reg min_temp post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72
matrix define beta = e(b)
matrix define beta_8 = beta[1,1..3]

quiet: reg min_temp post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72
matrix define beta = e(b)
matrix define beta_9 = beta[1,1..2]

quiet: xtreg min_temp post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72, fe
matrix define beta = e(b)
matrix define beta_10 = beta[1,1..3]

quiet: xtreg min_temp post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72, fe
matrix define beta = e(b)
matrix define beta_11 = beta[1,1..3]

quiet: xtreg min_temp post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=52 & min_temp <=72, fe
matrix define beta = e(b)
matrix define beta_12 = beta[1,1..2]

matrix define coef_`s' = (beta_1, beta_2, beta_3, beta_4, beta_5, beta_6, beta_7, beta_8, beta_9, beta_10, beta_11, beta_12)

}

matrix define beta = (coef_1)

forvalues s = 2/10000 {
	matrix define beta = (beta \ coef_`s')
}

clear
svmat beta
matrix drop _all

rename beta1 cons_change_1
rename beta2 cons_pre_trend_1
rename beta3 cons_post_trend_1
rename beta4 cons_change_2
rename beta5 cons_pre_trend_2
rename beta6 cons_post_trend_2
rename beta7 cons_change_3
rename beta8 cons_trend_3
rename beta9 cons_change_4
rename beta10 cons_pre_trend_4
rename beta11 cons_post_trend_4
rename beta12 cons_change_5
rename beta13 cons_pre_trend_5
rename beta14 cons_post_trend_5
rename beta15 cons_change_6
rename beta16 cons_trend_6
rename beta17 min_change_1
rename beta18 min_pre_trend_1
rename beta19 min_post_trend_1
rename beta20 min_change_2
rename beta21 min_pre_trend_2
rename beta22 min_post_trend_2
rename beta23 min_change_3
rename beta24 min_trend_3
rename beta25 min_change_4
rename beta26 min_pre_trend_4
rename beta27 min_post_trend_4
rename beta28 min_change_5
rename beta29 min_pre_trend_5
rename beta30 min_post_trend_5
rename beta31 min_change_6
rename beta32 min_trend_6

save Temp\Bootstrap_Point_Estimates.dta, replace

** The standard deviations for the following point estimates represent the
** bootstrapped standard errors presented in columns 1-6 of Table 2.

sum cons_change_1 cons_pre_trend_1 cons_post_trend_1
sum cons_change_2 cons_pre_trend_2 cons_post_trend_2
sum cons_change_3 cons_trend_3
sum cons_change_4 cons_pre_trend_4 cons_post_trend_4
sum cons_change_5 cons_pre_trend_5 cons_post_trend_5
sum cons_change_6 cons_trend_6

** The standard deviations for the following point estimates represent the
** bootstrapped standard errors presented in columns 1-6 of Table A9.

sum min_change_1 min_pre_trend_1 min_post_trend_1
sum min_change_2 min_pre_trend_2 min_post_trend_2
sum min_change_3 min_trend_3
sum min_change_4 min_pre_trend_4 min_post_trend_4
sum min_change_5 min_pre_trend_5 min_post_trend_5
sum min_change_6 min_trend_6
